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Context. The cosmic microwave background (CMB) spectrum probes physical processes and astrophysical phenomena occurring at 
various epochs of the Universe evolution. Current and future CMB absolute temperature experiments are aimed to the discovery of the 
very small distortions such those associated to the cosmological reionization process or that could be generated by different kinds of 
earlier processes. The interpretation of future data calls for a continuous improvement in the theoretical modeling of CMB spectrum. 
Aims. In this work we describe the fundamental approach and, in particular, the update to recent NAG versions of a numerical code, 
KYPRIX, specifically written for the solution of the Kompaneets equation in cosmological context, first implemented in the years 
1989-1991, aimed at the very accurate computation of the CMB spectral distortions under quite general assumptions. 
Methods. We describe the structure and the main subdivisions of the code and discuss the most relevant aspects of its technical 
implementation. 

Results. We present some of fundamental tests we carried out to verify the accuracy, reliability, and performance of the code. 
Conclusions. All the tests done demonstrates the reliability and versatility of the new code version and its very good accuracy and 
applicability to the scientific analysis of current CMB spectrum data and of much more precise measurements that will be available in 
the future. The recipes and tests described in this work can be also useful to implement accurate numerical codes for other scientific 
purposes using the same or similar numerical libraries or to verify the validity of dilferent codes aimed at the same or similar problems. 
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1. Introduction 

The CMB spectrum emerges from the thermalization redshift, 
Ztherm ~ 10^-10^, with a shapc very close to a Planckian 
one, owing to the tight coupling between radiation and matter 
through Compton scattering and photon production/absorption 
processes, radiative Compton and bremsstrahlung. These pro- 
cesses were extremely efficient at early times and able to re- 
establish a blackbody (BB) spectrum from a perturbed one 
on timescales much shor ter than the expansion time (see e.g. 
Danese & de Zottil (Il977b ). The value of Zthem dBurigana et alj 
1991 a') depends on the baryon density parameter, Q^, and the 
Hubble constant, Hq, through the product D - Q.b{H()/5Q)^ (Hq 
expressed in Km/s/Mpc). 

On the other hand, physical processes occurring at red- 
shifts z < Ztherm may leave imprints on the CMB spectrum. 
Therefore, the CMB spectrum carries crucial informations on 
phy sical processes occur ring during early cosmic epochs (see 
e.g. iDanese & Burigaria] (|T993) and references therein) and 
the comparison between models of CMB spectral distortions 
and CMB absolute temperature measures can constrain the 
physical pararneters o f the considered dissipation processes 
(lBuriganaetanil991al) . 
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The timescale for the achievement of kinetic equilibrium be- 
tween radiation and matter (i.e. the relaxation time for the photon 
spectrum), tc, is 



tc = t. 



ye 



4.5 X 1Q^\Tq/2.7 Ky^^-^Q-\l + zT'^ sec, 



(1) 



where fy,. - l/{neO-Tc) is the photon-electron collision time, 
tp = (Tg/Tr), Tg and = ro(l + z) being respectively the elec- 
tron and the CMB radiation temperature; kTe/nieC^ (being rrie the 
electron mass) is the mean fractional change of photon energy in 
a scattering of cool photons off hot electrons, i.e. » T/, To is 
the present radiation temperature related to the present radiation 
energy density by CrO - slTq (here a - Snlik*/ (hcf, h - jr"*/ 15). 
A primordial helium abundance of 25% in mass is assumed in 
these numerical estimates. It is useful to introduce the dimen- 
sionless time variable jeiz) defined by 



yeiz) 



r'" dt _ ri+' 

Jl tc ~ Ji 



exp 



dil+z) t. 
l+z tc 



(2) 



where t is the time, fo is the present time, and texp = l/H = 
l/[{da/dt)/a] is the expansion time, a = 1/(1 + z) being the 
cosmic scale factor normalized to the present time. 

As particular cases, by neglecting the cosmological constant 
(or dark energy) contribution we have 



6.3 X 10 
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where Zeq = 1.0 x \Q\TqI2.1 K)-^^„y is the redshift of equal 
non relativistic matter and photon energy densities and k - 
1 + A^v(7/8)(4/ll)'*/^ Ny being the number of relativistic, 2- 
component, neutrino species (for 3 species of massless neutri- 
nos, K - 1 .68), takes into account the contribution of relativis- 
tic neutrinos to the dynamics of the Universe |, while assuming 
Qa: - 0, Oa = 1 - Q.„r, and neglecting the radiation energy 
density, as possible at relatively low redshifts, we have 



-1/2 



sec , 



(4) 



where I/Hq - 3.1 x lO'^/i"' sec (h = //o/lOO). 

The time evolution of the photon occupation number, 
T]{v, f), under the combined effect of Compton scattering and 
of photon production processes, namely radiative Compton 
(RC) fcould 1984), bremsstrahlung (B) (Karzas & LattiB 
Il96ll: [Rybicki & Ligh tmanI Il979t) . plus other possible pho- 
ton emission/absorption contributions (EM )^ is well described 
by the complete Kompaneets equation (iKompaneetsI Il956t 
iBurigana et al|[T995l) : 



di 



1115 

<p tc dx 



+71(1 +7]) 

OX 



(5) 



dr] 




drj 




drj 




dt 


+ 

RC 


'dt 


+ 

B 


'di 


EM 



This equation is coupled to the time differential equation govern- 
ing the electron temperature evolution for an arbitrary radiation 
spectrum in the presence of Compton scattering, energy losses 
due to radiative Compton and bremsstrahlung, adiabatic cooling, 
plus possible external heating sources, q - a^^(dQ/dt), 



die 
dt 



(27/28)f, 



ey 



2T, 



dT, 



dt 



(32/27)^ 
3nek 



(6) 



here Teq,c = [hjr](l+ T])v'^dv]/[4k f yvVv] i s the 
Compton equilibrium electron temperature (iPeyraudI Il968t 
[Zeldovich & Levich 19_7Q), t^y - 3meC jAa-jer, - e^oCl + zf 
being e^o — ciTq the radiation energy density today, and x - 
hvo/kTo = /ivo(l + z)lkT(){\ -H z) is a dimensionless, redshift in- 
dependent, frequency (vq being the present frequency). 



2. Setting up the problem 

Partial differential linear equations are divided in three classes: 
elliptic, parabolic and hyperbolic. The Kompaneets e quation is a 
parabolic partial differential equation (lTricomilll957l) . Solutions 
to this equation under general conditions have to be searched 



'Strictly speaking the present ratio of neutrino to photon energy 
densities, and hence the value of k, is itself a function of the amount of 
energy dissipated. The effect, however, is never very important and is 
negligible for very small distortions. 

process that in principle should be included is the cyclotron 
emission. On the other hand, for realistic values of cosmic magnetic 
field, the cyclotron process never plays an important role for (global) 
CMB spectral distortions when ordinary and stimulated emission and 
absorption are properly taken into account and C MB realistic distorted 
spectra are considered (Zizzo & Buriganj '2005). In fact, the cyclotron 
term may be significant, in the case of deviations of ;; from the BB 
distribution at the electron temperature, only at very long wavelengths, 
corresponding to the cyclotron frequency, where, during the formation 
of a spectral distortion, FF and RC are able to keep rj extremely close to 
the BB equilibrium. 



numerically, because it is impossible to find analytical solutions 
that accurately take into account the many kinds of cosmological 
scenarios and the great number of relevant physical processes. 
The numerical code KYPRIX was written to overcome the lim- 
ited applicability of analytical solutions and to get a precise com- 
putation of the evolution of the photon distribution function for 
a wide range of cosmic epochs and for many cases of cosmolog- 
ical interest (Burigana et al. 1991a). KYPRIX makes use of the 
NAG libraries dNAG Ltdll2009l) . 

Besides these libraries, a lot of numerical algorithms are used 
in the code: we used some of the routines available to the scien- 
tific community, but often we did write routines dedicated to a 
specific task. Among the formers, the D03PCF routine of the 
current version of the NAG release has been used to reduce the 
Kompaneets equation into a system of ordi nary differential equa- 
tions (Dew & Wals'ij ll981l: iBerzins et al.| fl989: Berzins 1991 
Skeel & Berzins 1990). The D03PGF routine used in the first 
versions of KYPRIX is no longer available (see also Sect. l3.2l i. 
The two codes work adopting the same numerical framework or, 
in other words, the D03PCF routine of the current NAG release 
corresponds to the D03PGF routine of the NAG release used in 
the first versions of KYPRIX. On the other hand, they come from 
different technical implementations and exibit remarkable differ- 
ences in several aspects. The main differences between the two 
routines, their usage, and the corresponding implications for the 
code KYPRIX will be described in this work. In order to use the 
D03PCF routine, we have to put the Kompaneets equation in the 
form 



NPDE 



7=1 



(7) 



where Y is the time variable and X the spatial variable. Here 
/ = 1 and NPDE = 1. P^, Qi,Ri depend on x, t, U, dU/dx, the 
vector U is defined as immediately below. In our case, Pij - 1 
and m = (Cartesian coordinates). Note that Pij, Qi,Ri do not 
depend on dU/dt. 

The variables that enter in this equation are introduced and 
used in logarithmic form, to have a good and essentially uni- 
form accuracy of the solution in the whole considered frequency 
range. They are X-log(x) and U-log{T]). Hereafter, we will use 
the variables X, U, Y for what concerns the informatic aspect of 
the problem, keeping the use of the variables x, tj, y for consider- 
ations directly linked to physical aspects. The function Ri is de- 
termined only by the inverse Compton term while the other phys- 
ical processes, i.e. at least Compton scattering, bremsstrahlung, 
and radiative Compton, are included in the function Qi. In order 
to reduce Eq. ^ into a system of ordinary differential equations, 
the D03PCF routine uses the method of lines: the right mem- 
ber of Eq. (Q is discretized, reducing the calculation of partial 
derivatives in terms of finite values of the solution vector U at 
all the points of the X axis grid. Sp atial discretization is mad e 
by the method of finite differences dMitchell & Griffithsiri980l) . 
The resulting system of ordinary differential equations is solved 
using a differentiation formula method. The choice of the time 
parameter was driven by the need to have a very simple form of 
the Kompaneets equation. Finally, a "temperature independent" 
(time) Comptonization parameter 



Y = y(t) 



CTjC Trdt , 



(8) 



has fo und to be particularly advantageous dBurigana et al.l 
fT99Tah . 
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2.1. Boundary conditions 

Integrating equations of the type of Eq. (|7]) means to calculate 
the time evolution of the function U (X, Y), for a given initial 
condition U {X, 0) (in fact, the problem is also called "problem 
at initial conditions"). Numerically, the derivatives of U are re- 
placed by finite differences between values of U computed at a 
grid of points (in X) and the differential equation is replaced by 
a system of more simple equations. However, in presence of the 
only initial condition, this system is singular (Press et al. 1992). 
For this reason, resolving partial diff'erential parabolic equations 
needs boundary conditions: the problem is at initial values for 
the Y variable and at the boundary values for the X variable. In 
general, boundary conditions mean additional relations written 
to be joined to the system derived from the discretization to fi- 
nite diff'erences. 

Therefore, a good statement of the problem needs the defini- 
tion of appropriate boundary conditions. The capability of a 
refresh of these conditions along the integration in time leads 
more stability to the solution evolution because of the evolu- 
tion of the radiation field. Thanks to the opportunity of hav- 
ing the correct value of (p for each time step, the update of 
the boundary conditions can be physically motivated (see also 
Sect. |3.2.6t . The limits of the frequency range considered are: 
Xmin = log(jc„„„) = -4.3 e X,„ax = log(x„,a^) = 1.7. Of coursc, 
we want a solution of the Kompaneets equation over all the fre- 
quency range where it is possible to measure the CMB. Also, we 
need to consider a frequency range large enough to contain, in 
practice, all the energy density of the cosmic radiation field. 

The frequency range is so wide for the other two reasons 
described below. 

During the time evolution, some spurious oscillations of the 
solution at points close to the boundaries may appear (these ef- 
fects, that could also occur independently of the need of refresh- 
ing (p - for example for cases at constant (p -, may be partially 
amplified if, for computational reasons discussed in the follow- 
ing, the necessary refresh of the electronic temperature is not 
made for every time step). Fixing the frequency integration range 
limits far from the interval where we are interested to compute 
the photon distribution function allows to prevent the "contam- 
ination" of the solution by this possible spurious oscillations in 
the frequency range of interest. 

Since we can generally assume that a Planckian spectrum at 
Xmin is formed before recombination in a timescale shorter than 
the expansion time and, on the contrary, at x,„„^ the shape of 
the spectrum is unknown, it has been implemented in the code 
the possibility to adopt a particular case of Neumann boundary 
conditions: the requirement that the current density, in the fre- 
quency space, is null at the boundaries of the integration range 
(IChang&Coopeiill970l) : 



dx 



+ 77(1 +rf) 



(9) 



This choice of boundary conditions formally satisfies the re- 
quirement of the problem when we integrate the Kompaneets 
equation in the case of Bose-Einstein like distortions (with a 
frequency dependent chemical potential, fi - fiix), vanishing at 
very low frequencies). In fact, such distorted spectra are indis- 
tinguishable from a blackbody spectrum at sufficiently high and 
at low frequencies. 

Of course, it is possible to make a different choice of the 
boundary conditions by selecting Dirichlet like conditions. In 
this case the photon occupation number at the boundaries of 



the integration interval does not change for the whole integra- 
tion time. (In general cases, keeping constant conditions at the 
boundaries could be dangerous for the continuity of the solu- 
tion. Nevertheless, for some specific problems this condition can 
work - typically for problems with constant (p). 



3. A detailed view on KYPRIX 

The code KYPRIX has been written to solve the Kompaneets 
equation in many kinds of situations. The physical processes 
that can be considered in KYPRIX are: Compton scattering, 
bremsstrahlung, radiative Compton scattering, sources of pho- 
tons, energy injections without photon production, energy ex- 
changes (heating or cooling processes) associated to <^ 1 
at low reds hifts, radiative decays of m assive particles, and so 
on (see e.g.'Panese & Burig ana I (11993 ') for some applications). 
This code could be easily implemented to consider other kinds 
of physical processes. Various kinds of initial conditions for 
the problem can be considered and many of them have been 
akeady implemented in KYPRIX. The first obvious case is a 
pure Planckian spectrum. Several ways to model an instanta- 
neuos heating implying deviations from the Planckian spectrum 
have been introduced: a pure Bose-Einstein (BE) spectrum or 
a BE spectrum modified to become Planckian at low frequen- 
cies (this option could be exploited to integrate the Kompaneets 
equation with a constant (p and constant boundary conditions); a 
grey -body spectrum; a superposition of blackbodies. 

The data are saved into five files. 
DATI. This file contains the information about the specific 
parameters of the considered problem with a general description 
of its main aspects. 

DATIP. In this file we give the evolution of interesting quantities, 
like time, redshift, (p, and another quantities inherent to physical 
and numerical aspects of the problem (see also Sect. l5.lt . 
DATIG. It contains: the grid of points for the X axis used by 
the main program (remember that we are using a dimensionless 
frequency), a Planckian spectrum at temperature Tq and the 
solution vector U (that is to say log(?/)) at y = (starting time). 
DATIDE. This is the fundamental output file: it gives the 
solution of the Kompaneets equations at the desired cosmic 
epochs. 

DATIT. It is similar to the file DATIDE, but it contains the 
solution in term of brightness temperature (i.e. equivalent 
thermodynamic temperature; see Eq. ( fT2l i in Sect. l3.2.2] l. 



3.1. Main subdivisions 

The code is divided in several sections and, from a general point 
of view, is structured as described here below. 

1. Main program, in which many actions can be carried out: 
choice of the physical processes, choice of the cosmological 
parameters, initial conditions, characteristics of the numerical 
integration (accuracy, number of points of the grid), time 
interval of interest, choice of the boundary conditions, chemical 
abundances, ionization history. 

2. Subroutine PDEDEF. It is the subprogram where the problem 
is numerically defined. This subroutine is also divided in 
subsections to allow modifications in a simple and practical 
way. 

3. Subroutine BNDARY. Here the boundary conditions are 
numerically specified. 

4. Subroutines and auxiliary functions to perform specific 
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calculations. 



3.2. Technical specifications and code implementation 

The first version ^ of KYPRIX worked with the Mark 8 version 
of the NAG numerical library and were based on the routine 
D03PGF. The version of the NAG numerical library currently 
distributed is the Mark 21. Therefore an update of the code 
KYPRIX is necessary to adapt it to this new package. 

When KYPRIX starts running it asks all the input data: from 
the declarations of the output files' names to the integration ac- 
curacy and features. In the following subsections we give a de- 
scription of the various aspects of the code (and of its update), 
trying to give relevant hints about computational aspect of the 
code. 



3.2.1. Grid 

The frequency integration interval is divided in a grid of points 
(the mesh points): larger the number of points smaller the 
adopted frequency step. We adopted an equispaced grid in X. 
It is possible to used a very dense grid (for example 36001 mesh 
points corresponding to 36000 frequency steps). In general, it 
is necessary to use at least 3001 mesh points to have a solution 
accurate enough. 

We found an important difference between the two NAG ver- 
sions, not reported in the documentation of the routine D03PCF. 
In the first version (D03PGF), the subroutine where the partial 
differential equation is defined adopted the same mesh points de- 
fined in the main program. In the Mark 2 1 version the calculation 
is carried out in a different manner: the mesh points used in the 
subroutine PDEDEF is shifted of half spatial step with respect 
to the mesh defined in the main program. In this way, the mesh 
points in the PDEDEF subroutine will be exactly in the middle 
of the steps defined in the grid of the main program. For this rea- 
son, the limit of the integration interval are not considered in the 
mesh points in the subroutine PDEDEF and they are used only 
for the boundary conditions. 

The effect of this feature implies the definition of new parame- 
ters that play a fundamental role in the subroutine PDEDEF. The 
integral quantities in the Kompaneets equation (necessary to de- 
fine the radiative Compton term in the kinetic equations and the 
electron temperature) are computed once for any time step, in- 
side the PDEDEF subroutine. For this computation, arrays of 
dimension equal to the number of mesh points of the x variable 
as defined in the PDEDEF subroutine are used. Therefore, a par- 
ticular care must be taken in the definition of the dimension of 
the arrays defined in KYPRIX. Those used in the main program 
have dimension equal to the number of points of the mesh de- 
fined in the main program. The same dimension is given for the 
arrays defined for the boundary conditions. On the other hand, 
the major number of arrays are used in the PDEDEF subroutine 
to compute the integral quantities. The "inner" grid adopted in 
the PDEDEF subroutine is based on mesh points in the middle of 
the spatial steps of the main program grid, so the two grids can 
not work with the same point number; in fact, the arrays used in 
the PDEDEF subroutine have dimension NPTS - 1 . Therefore, 
in the main program and in the subroutine BNDARY we have to 
work with arrays based on the formula: 

X(I) = A H- (/ - 1) X , (10) 

V ; V ; (NPTS - 1) 



with I <X< NPTS , 

to define the correspondence between the grid of NPTS points 
and the X position, while we need another expression able to 
shift of half step the grid in the PDEDEF subroutine and based 
on NPTS - 1 mesh points: 



A + (I-\)x 



(B-A) 
(NPTS - 1) 



(B-A) 
liNPTS - 1) j 



(11) 



with \<X< NPTS - 1. 



For continuity reasons, we need to define (according to the 
choices made in the main program) the solution vector, contain- 
ing the photon initial distribution function, at the beginning of 
the integration also according to this grid definition. This vector 
is used by the PDEDEF subroutine as initial spectrum adopted 
for the computation of the rates of the physical processes and, 
of course, it is then renewed at every time step incrementation. 



3.2.2. Output 

Concerning the output files, the update version of KYPRIX 
stores a new vector containing the "inner" X grid used by the 
PDEDEF subroutine, XXGR (XGR refers to the main program 
Xgrid). 

In addition, we preferred to have the possibility to perform 
the conversion of the solution into equivalent thermodynamic 
temperature directly into the code and save it in a new output file 
(DATIT). The conversion relation is: 



xTq 



- term.eqinv 



ln(l + 1//7) 



(12) 



^Written in 1989 by C. Burigana. 



(we remember that in the code X - log[g(ji:) and U - logjo(7;)). 
The fundamental reason to perform this conversion directly in 
the code is associated to the extreme accuracy required for the 
solution in the case of very sm all distortions, of p articular in- 
terest given the FIRAS results (Fixse n et al.lll996l) . During the 
first tests, the conversion of the solution in brightness tempera- 
ture was performed at the same time of the solution visualiza- 
tion, through the IDL visualization program. The saving of the 
solution into files is typically performed not for all the points of 
the grid but for a reduced grid of, for example, 300 equidistant 
points along the original grid to avoid to store files of large size, 
useless for our scope, given the interest for the CMB continuous 
spectrum (by definition, the Kompaneets equation is not appro- 
priate to treat recombination lines). If the considered distortions 
were very small then the solution at each specific "inner" grid 
point could be affected by a numerical uncertainty not negligible 
in comparison with the very small deviations from a Planckian 
spectrum relevant in this cases. This numerical error is greatly 
reduced (becoming negligible for our purposes) by the averag- 
ing over a suitable number of grid points. Of course, the storing 
of the solution directly on a limited number of grid points makes 
this averaging no longer possible on the stored data. It were 
then necessary to average the solution values in intervals cor- 
responding to the output X grid directly into the code. Anyway, 
in many circumstances the diagram shape derived applying the 
conversion to brightness temperature only on the stored averaged 
solution still deviates at high frequencies from the effectively 
computed solution displayed by considering all the "inner" grid 
points because of the high gradients in the photon distribution 
function and/or in the brightness temperature that makes diffi- 
cult, or impossible, to find a general rule for the solution binning 
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that simultaneously works properly for the two solution repre- 
sentations. This problem is avoided converting the solution vec- 
tor in equivalent thermodynamic temperature before of the bin- 
ning of its values and then applying the binning to the equiva- 
lent thermodynamic temperature. The result is then a brightness 
temperature diagram very clean and precise, even for very small 
distortions. 

Other minor changes are made about the output data, where 
we passed from real to double precision, and for the saving fre- 
quency into the output files. 

3.2.3. Equation formalism 

A necessary update of the code has been performed to adapt it to 
the different formalism adopted by the new version of the NAG 
routine. This regards the expression of the Kompaneets equation 
in the PDEDEF subroutine. In particular, the D03PGF routine 
adopted the following expression of the partial differential equa- 
tion: 



and 



' dY ^ dX 



dX 



(13) 



where / - 1,2, ...,NPDE (number of partial differential equa- 
tions); C,, Fi depends on X, Y, U, dU/dX; Gtj depends on X, Y, U 
and U is the set of solutions values (Ui, U2, Unpde)- 

The expression now adopted by the D03PCF routine is rep- 
resented by Eq. (|7]l. 

It is simple to translate the code from the old to the new 
formalism. In the considered case NPDE = 1 . In this case, we 
have simply that R\ contains both the function Gi and the vector 
solution derivative with respect to X according to: 



7?! =Gii X 



dUi 
'dX 



(14) 



At this point, it is necessary to apply only the following sub- 
stitutions: 



Qi = -^'i 



and 



(15) 



With respect to this formalism, it is not difficult to 
adapt the various terms of the Kompaneets equation to the 
D03PCF routine. The terms that describe Radiative Compton, 
bremsstrahlung, (optional) electromagnetic processes and part 
of the contribution of the Compton scattering are counted in the 
function Q\. Instead, the second derivative of the solution vector 
with respect to X, which represents part of the inverse Compton 
rate, is counted in the function R\. So, according to these set- 
tings, we can write: 



Qx^FC + FBREM + FRAD + FDEC . 



(16) 



where FC stands for the contribution of the Compton scattering, 
FBREM the bremsstrahlung one, FRAD the radiative Compton 
one and FDEC represents the contribution of radiative decaying 
of particles. These terms are written in the form: 



FC = 



dx\dx / \dX 



-H 2 X 10^10^1 — -H 2 
dX 



dU 



In(lO) 



(17) 



FBREM + FRAD 



1 



X iFFO xWx 1.50"'^^ X FGAUNT + 



DCO 



W 



■hxGDC^ , (18) 



where FFO and DCO are the coefficients for the rates 
of bremsstrahlung and radiative Compton respectively; 
FGAUNT and GDC represent the Gaunt factor corrections 
for bremsstrahlung and radiative Compton, respectively; /i 
is an integral quantity refreshed at every time step (see also 
Sect. [323). 

The inverse Compton contributes to Ri in this way: 



dU 

/?! = — x0xln(lO) 
dX ^ 



-2 



(19) 



Finally, we can put Pi 1 - 1. 

Furthermore, it is possible to add other source terms'! in 
Eq. (O. 



3.2.4. Boundary conditions 

Also notable are the differences between the input expressions 
defining the boundary conditions. The D03PGF routine adopted 
an expression of the form: 

Pi(Y)Ui + Qi(Y)^ = Ri(Y, U) , (20) 

where / = 1 , 2, NPDE and Pi{Y), Rff, U), 2,(7) are functions 
to be defined. A quite different notation is used to provide the 
boundary conditions in the D03PCF routine: 



MX, Y)Ri{X, Y, U, Ux) = 7i(X, Y, U, Ux) , 



(21) 



where / = 1, 2,..., NPDE and /3i{X,Y)Ri(X,Y, U, Ux) and 
jiiX, Y, U, Ux) are functions to be defined (Ux = dU/dX). 
As a consequence of this notation, Neumann like boundary con- 
ditions can be now specified according to the expression: 



m = 1 

yd) = -XVA X (10^^'^ -H 1) X lnl0"2 



(22) 
(23) 



where XVA is the vector related to the X position computed in 
A and the dimension of both the equations corresponds to the 
differential equation number. Similar conditions are defined for 
the other extreme of the integration interval [A,B]. 



3.2.5. Accuracy parameters 

Another considerable difference between the two library ver- 
sions regards the definition of the integration accuracy param- 
eter D03PGF used three parameters for monitoring the local er- 
ror estimate in the time direction, supplying a good versatility. 
RELERR and ABSERR were respectively the quantity for the 
relative and absolute component to be used in the error test. The 
third parameter, INORM, was used to define the error test. If 
E(/, j) is the estimated error for t/,- (the vector solution) at the 
j - th point of the X grid, then the error test was: 



''For example, an already implemented subroutine models a term, 
called here FDEC, to be added in Eq. l ll6b which acconts for contribu- 
tions of possible radiative decays of massive particles in the primordial 
Universe. 
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- INORM = ^1 E(/, |< ABSERR + RELERRx | U{i, j) \ 

- INORM = 1 ^1 E(/, j) |< ABSERR + RELERR x max,. | 
U{i,j)\ 

- INORM = 2 ^ ||E(i,;)|| < ABSERR + RELERR x ||f/(/,;)||. 

Instead, according to the new library version we have to define 
only one parameter ACC, a positive quantity that monitors the 
local error in the time integration. If E(/, j) is defined as above, 
then the error test is: 

|E(/,;)|=ACCx(l+|f/(/,;)|). (24) 

Note that this is equivalent to the error test implemented in 
the D03PGF routine in the case INORM = and ABSERR = 
RELERR (= ACC). 

3.2.6. Electron temperature 

During the numerical integration, some subprograms use the dis- 
tribution function calculated at that time to compute 4>- The inte- 
grals to be computed are those that we find in the expression for 

'Peq,C- 

'Peq,c = = — TTSS — ■ (25) 

In this calculation, the integration range is obviously the inte- 
gration interval considered for the problem: A < X < B (that, 
in terms of mesh ordering, corresponds to the range between 1 
and NPTS or NPTS - 1). For computing these integrals, all the 
points of the grid are used. The integration is based on the NAG 
DOIGAF routine, suitable for tabulated functions. The update 
value of (p is also used in the boundary conditions. 

In the previous version of the code KYPRIX, the compu- 
tation of integral quantities were performed through a specific 
modification of the NAG package implemented by the KYPRIX 
code author that allowed to recover the whole vector solution at 
each time step in the subroutines (and in particular in PDEDEF), 
while the original package made only available in PDEDEF the 
solution separately at each grid point (being in fact the pack- 
age originally designed for "pure" partial differential equation, 
without terms involving integrals of the solution). This modifi- 
cation, possible thanks to the availability of the NAG sources 
(and, in practice, thanks to the relative simplicity of the early 
library versions), permitted to update the integral quantities per- 
fectly according to the "implicit" scheme adopted by the code 
for the integration in time. This is no longer feasible. Therefore, 
the update of the integral quantities must be now performed with 
a "backward" scheme, saving the solution at the previous time 
step in a proper vector and using it in the computation at the 
given time step. As well known, "backward" schemes are typ- 
ically less stable than implicit schemes. And, in fact, we veri- 
fied in some cases the difficulty of the D03PCF routine to work 
implementing the update of the quantities coiTesponding to the 
integral terms in the Kompaneets equation (and in particular of 
cj)) for each time step. This were likely due to numerical instabil- 
ities. 

We have then introduced a new integer control parameter 
into the code: STEPFI. It determines the frequency for the update 
the dimensionless electron temperature 0, relevant, of course, in 
the case we want to perform an integration with a variable (p. We 
have checked that updating the integral terms in the Kompaneets 
equation not at every time step, but after a suitable number of 
time steps does not affect the accuracy of the solution. This is 
due to the fact that the time increasing in the code is performed 



with very small steps while the physical variation of occurs on 
longer timescales^. 

See Sect. 15.2.21 for tests regarding the implications of this 
new implementation of the electron temperature evolution. 

3.2.7. Radiative Compton 

In the computation of the radiative Compton term there is an 
integral term, /i , given by the numera t or of t he right member of 
Eq.|25](see Eq. (16) in lBurigana et"an dHH), so It IS necessary 
to harmonize its update according to the parameter STEPFI 
discussed in the previous subsection. In fact, a possible asyn- 
chronous update of it and (p could create numerical instabilities 
and the crash of the code run, as physically evident from the 
great relevance of both radiative Compton term (at least at high 
redshifts) and electron temperature for the evolution of the low 
frequency region of the spectrum. 



3.2.8. Integration routines 

The global accuracy of the code KYPRIX depends on the 
accuracy of the solver for the partial differential equation as well 
as on the accuracy of all the other routines dedicated to different 
specific computations. In this section we focus on the routines 
of numerical integration used in the KYPRIX. As discussed in 
the previous subsections the DOIGAF routine has been used 
for tabulated functions. On the other hand, KYPRIX involves 
the computation of integrals of various functions defined by 
analytical expressions, namely the parameters characterizing 
different types of initial conditions for the distorted spectra, 
as the amount of fractional injected energy and the electron 
temperature, and the relationship between the different time 
variables entering in the code (see also Sect. 14.11 ). Ultimately, 
the better accuracy in these specific computations implies a 
better global accuray of KYPRIX. In the early release of the 
code KYPRIX the NAG DOIBDF routine was used to calculate 
integrals of a function over a finite interval. The same task 
can be caiTied out by the DOIAJF routine. This code offers a 
better accuracy than DOIBDF (DOIAJF is in fact suitable also 
to integrate functions with singularities, both algebraic and 
logarithmic). After the routines substitution, the results showed 
a great increasing of accuracy. In particular, this improvement 
offers the possibility to investigate also very small distortions 
that requires a very precise determination of all the relevant 
quantities because the absolute numerical eiTor of the integration 
must be much smaller than the (very small quantities) of interest 
in these cases. In particular, the quantity Ae^/e, (where is the 
actual density energy and e, is the energy density corresponding 
to the unperturbed distribution function just before the energy 
injection) must be constant during all the integration process in 
the absence of energy injection terms, according to the energy 
conservation. The precision increase on the computation of 
this quantity was noteworthy, keeping now always inside a 
few percent of the physical value (and its possible physical 
variation; see Sect. 15. It of the same quantity independently 
of the magnitude of the considered distortion, allowing at the 
same time to accurately check the global accuracy of the code, 
improved thanks to the better computation of all the integral 



'Of course, for physical processes with a stronger variation of the 
electron temperature, the accuracy parameter (see previous subsection) 
should be good enough to force the code to adopt sufficiently small time 
steps. 
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terms appearing in the Kompaneets equation. The remarkable 
improving in energy conservation even in the case of very 
small distortions is an important feature of the new version of 
KYPRIX that makes is applicable to a wider set of cases. 



properly evaluated considering also possible energy injections 
after neutrino decoupling). 

After some calculations, we can finally write the completed 
and updated expression of df/dw: 



4. New physical options 

4.1. Cosmological constant 

About ten years ago, the relevance of the cosmological con- 
stant term (or of dark energy contribution) has been probed by 
a wide set of astronomical observations of type la Supernovae 
dPerlmutter et al.lil999l:lRiess et al.lll998h . We have then updated 
the numerical integration code KYPRIX to include the cos- 
mological constant in the terms controlling the general expan- 
sion of the Universe. In particular the input background cosmo- 
logical parameters considered in the code are now: To,k,Ii[= 
i/o/(100Km/s/Mpc)],Q„,.,Qi,QA,QA:, i.e. the present CMB 
temperature, the contribution of massless neutrinos, the Hubble 
constant, the (non relativistic) matter and baryon energy density, 
the energy densities corresponding to cosmological constant and 
curvature terms. 

In order to compute the proper cosmic evolution of the var- 
ious term s, an (increasing with time) scale factor parameter co 
dsHE&Stebbins 1983), defined by 



a _ m^c 1 
fli kTo I + z 



1.98 X 10''©(1 +zy 



(26) 



has been adopted (here = To/3°K and the index 1 is referred 
to a particular epoch, when the CMB energy density was equal 
to the electron mass*]]; kjiai) - nigC^). To write a suitable ex- 
pression for its time evolution we have to introduce two new key 
parameters 



/? = ^=3.5xlO-''^Q„, 

Prl ®^ 



(27) 



that is the initial ratio between matter energy density and radia- 
tion energy density, and 



Tgl CO 
2.164x10" ■ 



2.164x10^' 



1/2 



(31) 



where Qj/,. = Q.x/Q.y. 

The equation for d> has to be inserted in the expression giving 
dy - fledf, which is inside the integral used to compute the time 
variable y(aj) = dy because we set y = when the integra- 
tion starts atoj - (Ostan (or equivalently at z = Zsian)- Finally, the 
expression for the time evolution of w and y are related by the 
variable change: 



dt o) 1 

dy = flcdf — fle — d(jj — Qc — — do) , 

d(jj <i) LO 



(32) 



where = (f>l{TciCjf) (r^i = 2.638 x 10"^ @^ l{h^Q.b)). 

The implementation of the cosmological constant and cur- 
vature terms makes the code KYPRIX suitable to be applied 
to interesting cases at late ages, including cosmic epochs at 
z < 1 when A supplies the greatest contribution to the expan- 
sion rate of the Universe. Remarkable examples are spectral dis- 
tortions associated to the reionization of the Universe that, in 
typical astrophysical scenarios, starts at relatively low redshifts 
(z i 10 - 20). For sake of completeness, a few words about the 

computation of the time evolution in the code. The subroutine 
called WDIYO is the core of the time evolution in KYPRIX: it 
computes the value of w given a value of y. This process takes 
advantages of the definition of y as integral of dy and, of course, 
it happens at each time step. To do this, we make use of a dou- 
ble precision version of the function ZBRENT, from "Numerical 
Recipes" (.Press et al. ,,1992) . 

4.2. Chemical abundances 



1 /Stt 

T„l \ 3 



1/2 



87r a 



„2\4 



1/2 



= 0.076s" 



(28) 



defined as a initial gravitational time scale. The quantities with 
the index 1 refer to the epoch when a - a\, with the index 
when f = f (today); p,-i and p„,i are related to pr and pm by 



In this new version of the code it is possible to choose the primor- 
dial abundances of H and He. We have consistently computed 
the electron number density, Wg, involved in the different physi- 
cal processes. In the previous implementation it was assumed a 
fixed abundance of H and He (and full ionization). The electron 
number density was then given by 



■■ POr\ 



1 

■Pr\—\ 

of 



Pm - PVm 



1 

■ Pm\ 



(29) 



free u 



Pb 7 

OTi 8 



(33) 



respectively. 

Now we can define an equation for the evolution of u: 



CO 
CD 



—Gp{(x)) 



1/2 



(30) 



Pr\K Pml PK\ 

or oj^ oj- 



■Ph 



where we have included the contribution of massless relativistic 
neutrinos in the term k (see also footnote 1 ; the term k should be 



*So, the parameter oi is analogous to the scale factor a, but normal- 
ized at the epoch in which a = ai, that is to say when kT = m^c^. 



where p/, is the baryon density and nib the mean mass of a 
baryon. Now, denoting with fn the fraction in mass of primordial 
H and considering that He - hh + 2nHe, we have 



free 



1 + /h Pb 
2 nil, 



(34) 



This obviously impacts the physical processes involved in the 
code. Note that Compton scattering, radiative Compton, and 
bremsstrahlung depend linearly^ on n{'^''^ . 



'For bremsstrahlung, the dependence on the densities of nuclei is 
now explicit. 
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4.3. Ionization history 

In the last years, CMB observations are getting a very high ac- 
curacy, in particular regarding the angular power spectrum of 
temperature anisotropics albeit also for that of E mode polariza- 
tion and cross-correlation (see e.g. iNolta et al] (|2009|) ). A de- 
tailed understanding of the cosmological reionization process 
is crucial for precisely modeling the power spectrum of CMB 
anisotropics in comparison with current data, while a better 
treatment of recombination is of great relevance in the view 
of the data expected by the forth coming ESA Planck satellite^ 
dThe Planck Collaborationll2006l) . 

An accurate modeling of the ionization history is crucial also 
for the precise computation of CMB spectral distortions in the 
view of the comparison with data from a future generation of 
high sensitivity experiments. We have then included these as- 
pects in the current implementation of KYPRIX. In particular, 
the fraction of each state of ionization of the relevant elements 
(hydrogen, H, and helium. He) has been implemented. 

The effects of this implementation is negligible^ in prac- 
tice for the radiative Compton, because during the epochs this 
process is active the medium is fully ionized. Instead, Compton 
scattering and bremsstrahlung rates are significantly influenced 
by this implementation. 

Once introduced the electron ionization fraction in the code, 
Xe, which gives the number of electrons that takes part in the 
physical processes, we can choose different ways by which the 
active fractions of elements can play their roles in the phenom- 
ena. 

Given Xe, from the charge conservation law we have a con- 
straint on the number of the free ions in the considered plasma. 
The simplest way to take count of them in the code is to assume 
an equal fraction of ionization for H and He. Of course, this is 
a toy model, but this parametrization was very useful to test the 
code. 

A somewhat more accurate treatment of the physics of reion- 
ization/recombination processes implemented in the code is 
based on the Saha equation: 



rii 



2 gi+i 



(35) 



where «, is the density of atoms in the i-th state of ionization, 
He is the electron density, g, is the degeneracy of states for the i- 
ions, e, is the energy required to remove / electrons from a neutral 
atom and A is the thermal de Broglie wavelength of an electron, 
defined by 



A = 



/z2 



^ ImngkBT 

Providing the electron ionization fraction defined as"]| 

„free 



Xe = 



(36) 



(37) 



* http://ww w.rssd.esa.int/planckl 

'For example, for z £ 10"'* and Ae,-/6, = 10"' we find that, for both 

BE and Comptonization like distortions, the rate of radiative Compton 
is less than 1/1000 of the rate of bremsstrahlung at each frequency of 
the grid. 

'"it is common to find in the literature normalized to the H num- 
ber density. Obviously, the code can switch between the two conven- 
tions. 



we can compute the ^xx^xvoyNxv'i,XH,XH*,Xne,XHe* ,XHe** , defined 
as the relative abundances of the different ionization states of 
each element with respect to the global number density of the 
element. The Saha equation provides the ratio between two state 
of ionization of a single specie, once given the electron density 
and the temperature: 



riHe 



nHe 



nue 



(38) 



To recover all the unknowns we need other relations. These 
additional conditions are provided by the charge conservation 
and the nuclei conservation. 

From charge conservation it is possible to recover the con- 
tribution of the electrons related to a single specie to the total 
number of free electrons. The latter is of course related to the 
ionization fraction of all the involved species and can be written 
as: 



Jree 



ph_ 

nib 



XH*fH + 2 I /" I (XHe* +XHe**) 



(39) 



where is the fraction of primordial H described in the previ- 
ous section. In Eq. (39[ it is possible to identify the number of 
electrons coming from H and from He. 

The nuclei conservation law separately for H and He is ex- 
pressed by 



'He 



Pb_ 

nit 

Pb_ 
nit 



[fHiXH +Xh*)] 



l-fH 



(XHe +XHe* + XHe^*) 



(40) 



(41) 



Providing the electron ionization fraction and the tempera- 
ture, it is possible to build two separate systems (one for H and 
one for He) in order to recover the considered ionization frac- 
tions. 

The code can ingest a table with the desired evolution of Xe 
and, as necessary for a physical modeling of cosmological reion- 
ization, of the electron temperature. 

Obviously, we have implemented the best way to perform the 
exact calculation of the rates of considered processes in scenar- 
ios involving reionization/recombination which is that of using a 
co-running code, coupled to KYPRIX, able to supply the ioniza- 
tion fraction for all the species. For the recombination process, 
we developed an interface that allows the code to call an external 
program, in our case RECFAST (Seager et al. 1999), and run it 
with the same cosmological parameters selected for KYPRIX. 
Then KYPRIX will read the output table from RECFAST and 
will use the ionization fractions for the calculation of the rates of 
the processes, allowing a more detailed estimation of the s pectral 
distortion arising during the recombination process ( Wong et al.l 
l2008h . If we are interested in standard recombination, we adopt 
the equilibrium electron temperature (Eq.lZSll. 

A consistent modeling of a modified recombination should 
provide the evolution of both electron temperature and ionization 
fractions (at least ;t'e). 

5. Tests on code performance, reliability, and 
accuracy 

Once terminated the updating of the numeric integration code, 
we have carried out many accuracy and performance tests. 
Obviously, we have checked the physical meaning of the nu- 
merical solutions provided by the code comparing them with the 
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existing analytical solutions that in some cases can be consid- 
ered as good approximations of exact solutions. In general, a 
code of good quality must have high numerical precision com- 
pared to the knowledge, both theoretical and observational, of 
the considered problem. Note also that various routines, mainly 
from the NAG pack age albeit also from the "Numerical Recipes" 
dPress et al. 1119921) package, are used here for several specific 
computations. Since different routines can solve the same math- 
ematical problems using different numerical methods and/or im- 
plementations with different settings and input parameters, we 
have verified that the adopted routines allow the appropriate ac- 
curacy and efficiency. 

In the following subsections we will report on some specific 
quahty tests of the numerical solutions. We comment now briefly 
on the code computational time. 

In order to evaluate the global CPU time of the code we per- 
formed many runs with very different settings. These times are 
carried out testing the code on a machine with 4 Alpha CPUs, 
but effectively using only one CPU (now we are running the code 
on IBM Powers Processors). The code global CPU time ranges 
from few minutes, for cases in which the integration starts at 
low redshifts, to about 5 hours, for some cases starting at very 
high redshifts {y(z) - 5). There are many variables that take role 
in determining the global CPU time. The complete Kompaneets 
equation involves in fact several terms. In the code KYPRIX we 
can select the physical processes to be considered in the numeri- 
cal integration and the global CPU time increases with the num- 
ber of activated processes. 

Of course, the global CPU time depends on the parameters 
related to the numerical integration characteristics. The number 
of points adopted for the X grid (see Sect. |3.2.TT l has a great influ- 
ence on the global CPU time. Clearly, the integration accuracy 
improves with NPTS and in many cases of interest it must be 
set to a very large value. We find that the global CPU time is 
approximately proportional to NPTS {tcpu °^ NPTS). 

The parameter that plays the most relevant role in determin- 
ing the global CPU time is the accuracy required for the time 
integration. The final solution precision depends on the value of 
the corresponding parameter ACC. Only for very high accuracy 
(ACC i 10 '2 - 10"*) the CPU time reaches the duration of 

some hours while keeping ACC ~ 10"^ the integration is carried 
out in few minutes. Anyway, the limits imposed by CMB spec- 
trum observations drive us to investigate in particular on small 
distortions. It is then necessary to work with low values of ACC 
(in general, ^ 10^^). 

After the assessment of the better choice for the use of the 
parameters of the various numerical routines and the definition 
of the accuracy of the integration in time, we have carried out 
several tests in order to verify the physical validity of the code 
results. 



5.1. Energy conservation 

In the code output file DATIP we store values of several param- 
eters of interest. Two of them provide very useful information 
on the goodness of the numerical integration. The first one is the 
ratio, er/ei, between the radiation energy density at each time 
and the energy density corresponding to the unperturbed distri- 
bution function before the distortion'! in the absence of dissipa- 



1.0000 



"For example, for a Bose-Einstein distorted spectrum with a di- 
mensionless chemical potential yUo produced by an energy dissipation 
with negligible photon production 6,-/e, = f(p.o)/>p'^^^(Mo) (see e.g. 




0.0001 



10* 10" 
redshift z 



Figure 1. Error in the energy conservation expressed in terms 
of relative (%) deviation from its input value of the quantity 
Aer/ei - 10"^. The accuracy of the integration in time is set 
to ACC= 10"'^ (see also the text for further details on com- 
putation). Solid, dot-dashed, and dashed lines refer to an en- 
ergy injection occurred respectively at ^^(z) ^ 1.5,0.25,0.01. 
Obviously, this error further decreases improving the accuracy 
parameter adopted in the numerical integration. 



tion processes, the exact energy conservation is represented by 
the Constance of this ratio during the integration. To quantify the 
accuracy supplied by KYPRIX, the values of er/f/ are stored at 
the starting of the integration and at many following times. In 
order to estimate the precision of the energy conservation, we 
define the quantity: 



ERR, = 



l(er/e,)/=/„„„ - (er/edt: 



1 



(42) 



that gives the relative induced error in the energy conservation 
with respect to the initial value of Ae^/e, or, more formally, the 
relative error induced by numerical uncertainty on the amount of 
fractional injected energy. Typical results obtained starting from 
a superposition of blackbodies are reported in Fig.[T] 

Since the same absolute numerical integration error corre- 
sponds to a larger relative error for a smaller distortion, i.e. for 
smaller Aer/e,- in these models, we could in principle expect 
a relative degradation of the energy conservation for decreas- 
ing distortion amplitudes. Our tests indicate in fact an increas- 
ing of the relative degradation of the energy conservation for 
decreasing distortions when the same accuracy parameters are 
adopted. On the other hand, one can select them according the 
specific problem. For example, considering an energy injection 
of Aer/ei - 10"^ occurred at z ~ few x 10"*, an accuracy equal 
to ACC = 10"^ is fully satisfactory for a very precise calculation 
of the photon distribution function, while for earlier processes, 
as forz > fewx 10^, and for the same injected fractional energy, 
the accuracy parameter needs to be set to ACC < 10"'^ in order 
to assure a calculation with ERR, less than 1%. Anyway, we find 
that for suitable choices of the integration accuracy parameters, 
ERR, can be kept always below ^ 0.05% without requiring a too 



iSunvaev & Zeld ovich ( 1970), Danese & de Zottil ( [19771) and also equa- 
tions (8) and (9) in Burigana et al.. (.1991a ) for the definition of /(yUo) 
and (/)(yio) for arbitrary values of //o; for //o <K l^/O^o) - 1-1.11/Joand 
V(fio) - 1 - 1.368;Uo). 
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large computational time. Finally, we note that in some circum- 
stances the scheme for the electron temperature evolution in the 
new version of KYPRIX (backward differences), different from 
that used in the original one (implicit scheme), could imply some 
small discontinuities in the evolution of the electron temperature 
and of Ae,./ e,. On the other hand, we have verified that this effect 
does not affect the very good accuracy of the solution, because of 
the very small amplitude of these discontinuities (together with 
their localization to a very limited number of time steps) and of 
the corresponding energy conservation violation. 

5.2. Comparative tests 
5.2.1. Comparing solutions 

The first kind of test is simple comparison between the results 
obtained with the update version of KYPRIX and those obtained 
with the original version for the same set of input parameters. 

To this purpose, we have considered some interesting cases 
carried out in the past. I n particular we used the input parameters 
adopted in jBurigana et al. 1995) where also semi-analytical de- 
scriptions of the numerical solutions of the Kompaneets equation 
were reported. We report here cases characterized by an amount 
of exchanged fractional energy Ae^-Ze, = 10""*. We started the 
integration from a redshift corresponding to yhiz) - 0.25 in one 
case and to yniz) - 0.01 for another. The input cosmological 
parameters are: Hq = 50, Q/, = 0.03, k = 1.68, To = 2.726 K. 
The results given by the update version o f KYPRIX (see Figj2]l 
are fully consistent with those reported in lBurigana et al.l(ll995h . 
Moreover, since that paper provided also a seminalytical de- 
scription of the solution of the Kompaneets equation, it is clear 
that a good agreement of the numerical results obtained with the 
original and update code represents a further confirmation of the 
validity of the analytical description. 



5.2.2. Electron temperature behavior 

While the energy exchange between matter and radiation is 
leaded by Compton scattering (in the absence of external energy 
dissipation processes), electrons reach the equilibrium Compton 
temperature (Pevraud 1968; Zeld ovich & LevichI 19701) . 7^ ^^, in 
a time shorter than the expansion time. For yi, «; 1, (late heat- 
ing processes) Compton scattering is no longer able to signif- 
icantly modify the shape of the perturbed spectrum and the fi- 
nal electronic temperature (pj (remember that = T^ITr), im- 
mediately after the decoupling, is very close to (p^q - T^ eg/Tr. 
On the other hand, for energy injections at redshifts correspond- 
ing to yi, > 5, Compton scattering can establish kinetic equilib- 
rium between matter and radiation. This corresponds to a Bose- 
Einstein spectrum, with a final electron temperature given by 
(Sunyaev & Zeldovich 1970; Danese & de Zottll977.) : 



0/(3'A>5) = 0ij£^(l-l.lljUo) 



-1/4 



(43) 



where //o(^ 1) is the usual dimensionless chemical potential. 
Moreover, in this case the evolution of the chemical potential, 
fj.(z), the relation between it and the amount of fractional energy 
injected, Ae/ e,, depends on the energy injection epoch. 

For the intermediate energy injection epochs, corresponding 
to yi, < 5, the final value of </> (a function depending on y/,) is 
between the values of (/)be and (p^q, because the Compton scat- 
tering works to produc e a Bose-Einstein like spectrum anyway 
dBurigana et alJll991ah . At these epochs the relation between 
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Figure 2. Comparison between the present time solution for the 
CMB spectrum obtained from the old version of the nu merical 
code (top panel; adapted from panel a of Fig. 1 in Burigan a et al.l 
(119951) ) and the current one (bottom panel). See the text for 
further details on computation parameters. Dashed (dot-dashed) 
line refers to an energy injection occurred at yh(z) - 0.01 
iyiiiz) - 0.25). Note the excellent agreement between the results 
of the two codes. In the case of the old version of the code we 
report also the analytical approximation (dotted line) described 
by a Comptonization spectrum plus a free-free distortion (see 
Eq.|50ll. 



the chemical potential and the amoun t of fractional injected en- 
ergy injected is simply given by CSunvaev & Zeldovichlll970l ; 
iDanese & de Zottilll977h : 



m ^ 1-4 ■ 



Ae 



(44) 



By exploiting the numerical results, a simple formula for <p can 
be found (Burigana et al. 1995): 



k 5 — yh 

<f>f(yh) = 7-; (Ac/ - 4>be) + <^B£ ; 

^ K + yh 



(45) 



here k - 0.146 is a constant derived from the fit. Moreover, this 
expression represents an accurate description of the evolution of 
for any value of yh- In facts, for a value of y (y < yh) they 
found: 



^iy,yh) = ^fiyii-y). 



(46) 



In the considered cases, as in many situations of interest, the 
perturbed spectrum of the radiation (immediately after the heat- 
ing process) is described by a superpositi on of blackbodies and 
the equilibrium temperature is given by dZeldovich & SunvaevI 
119691: IZeldovich et al.lfT972l:lBurigana et al.lll995l) : 



*eo - (1 +5.4m)0,-. 



(47) 



where 0,- = Ti/T,. = (1 H- Ae/e,)"'/"* ^ 1 - m and the 
Comptonization parameter u could be related to the am ount 
of fractional energy exchanged by (Zeldov ich & Sunyaey||T969l: 
IZeldovich et"aLlll972l:lBurigana et al...l995i) : 



M -(l/4)Ae/ei. 



(48) 
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Figures. Top panel: evolution of as derived from the numer- 
ical code. The parameters adopted for the computation are the 
same as in Fig. [T] as well as the adopted kinds of lines. We report 
also for comparison the analytical results obtained from Eq. ( l45T l 
(thin solid lines). Bottom panel: relative differences between the 
analytical results and the numerical results expressed in terms of 

{<!> analytical ' (pnmnerical) / (^numerical ' !)■ Note that in absolute Value 

they are < 1% at each time. Again, the adopted kinds of lines 
are as in Fig.[T] but solid line is replaced by dots when the above 
relative difference is negative. 



Eqs. ( |45] ) and ( |46b can be used to test the behavior of the 
electron temperature during the numerical integration of the 
Kompaneets equation carried out with the new code version. 
With the increasing of the time variable, the values of (p are saved 
into the file DATIP, from the initial time step to the final one. The 
upper panel of Fig. [3] shows the two behaviors of (f> (the numer- 
ical one and the analytical expression given by Eqs. ( |45] l and 
( |46] )) while the bottom panel displays their relative difference. 

This test is of particular importance for the check of the 
validity of the results because of the crucial role of (p in the 
Kompaneets equation. As remembered in the previous section, 
in the new version of the code a different scheme is used with 
respect to that implemented in the original version. The veri- 
fication of the very good agreement of the above behaviors of 
(p further supports the substantial equivalence of the two code 
versions and their reliability. In particular, it confirms that the 
new adopted numerical scheme for the evolution of 0, in princi- 
ple less stable than the implicit scheme implemented in the code 
original version, does not affect the validity of the solution. 

5.3. Free-free distortion 



As already discussed bv lSunvaev & Zeldovichl(ll970h . accurate 
measures of the CMB spectrum in the Rayleigh-Jeans region 
could provide quantitative informations about the thermal his- 
tory of the Universe at primordial cosmic epochs. On the other 
hand, photon production processes (mainly radiative Compton 
at earlier epochs and bremsstrahlung at later epochs) work to 
reduce the CMB sp ectrum depression at long wavelengths (see 
iDanese & de Zottil (l9m ) since they try to establish a true 



(Planckian) equilibriu m. For z < Zp dPanese & de Zottil 1 19801 : 
iBurigana et al.lll991ah with 



,,..,4XI0^(^)'"(4)'\.«, 



(49) 



low frequency photons are absorbed before Compton scattering 
moves them to higher frequencies. 

In case of small and late distortions (zi, < Zp), a quite good 
appro ximation of the whole spectrum is given by tBurigana et aTl 
[T995h 



-n) 



1 



dr' 



(gA-M _ 1)21 tanh(jc/20,) 



(50) 



where the index / denotes the initial value of the correspond- 
ing quantity and u is the Comptonization parameter This ex- 
pression provides also an exhaustive description of continuum 
spectral distortion generated in various scenarios of (standard 
or late) recombination or associated to the cosmological reion- 
ization. For an initial blackbody spectrum, at dimensionless fre- 
quencies xb x <g 1 the above equation can be simplified to 
jBurigana et al.ll 1 995l) 



^ yB 

VBBJ H T 



x/4>i 



(51) 



here yg, an optical depth of the Universe for bremsstrahlung ab- 
sorption (radiative Compton can be neglected at late epochs), is 
analogous to the Comptonization parameter and it is given by 



yB 



jti, 



((/>-4>i)'P'^^^8B{x,cf>)KoBdt 



(52) 



r 

Jl+Z 



/,-3/2 



gB(x, 4l)KQBtexp 



d(l + Z) 

l+z 



Xb is the frea i iency at which yabs,B - 1 dZeldovich et al.l 
il972t lde Zottil Il986l) The d ep endence of the Gau n t fac- 
tor ("Karzas & Lattei^ Il96ll: iRvbicki & LightmanI 119791: 
Burigana et al. 1991a) on x and (fi at very long wavelengths is 
weak: gB ln(x/0). 

In terms of brightness temperature, the distortions at low fre- 
quencies (at any redshift) can be written as 



Tbr - Tr4>i yB r, , 

- ^ - 2m<^; , 

1 r X'- 



(53) 



where T, = ro(l + z). This approximation holds at low frequen- 
cies but not at too low frequencies, where the brightness tem- 
perature obviously approaches the electron temperature because 
of the extremely high efficiency of bremsstrahlung, still able to 
generate a Planckian spectrum at electron temperature. 

In order to show that our numerical solution follows the be- 
havior described by the above equation, we can compute yB from 
the brightness temperature derived from the numerical solution: 



yB 



Thr - Tr4 
Tr 



+ 2u<pj 



(54) 



The reported numerical result (see Fig.|4]i refers to a heating pro- 
cess corresponding to a full reionization starting at z ^ 20 with 
<p = 10 X lO'* K, producing a final Comptonization parameter 
u - 4 X 10"'' compatible with FIRAS upper limits. As shown 
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in Fig. m at low frequencies ys approaches an almost constant 
value. This is correct only for A > 200 - 300 cm while at higher 
frequencies yg is no longer almost constant (see again Fig. |4|i 
because of the dependence of the Gaunt factor on x and <p as ex- 
pressed in the definition of ye- We can also write an expression 
describing the brightness temperature through a constant param- 
eter ys, derived from Eq. (|54l i (see Fig. |4]i averaged over the 
range at very long frequencies before the ys declining shown in 
Fig. in to verify the accuracy of the below expression 



Thr,ys = ( ^ - 2m</'/ + (pi] - T, 



(55) 



in comparison with the numerical results. Note that where the 
Gaunt factor dependence on x and produces a significantly 
varying ys (see Fig. HI, the Comptonization decrement is more 
relevant than the free-free excess in determining the brightness 
temperature, as evident from the good agreement of the two 
curves in Fig. H) Clearly, the brightness temperature derived in 
this way works only up to frequencies (~ 100 GHz) approaching 
that at which the excess in the brightness temperature produced 
by the Comptonization begins (at ~ 220 GHz, see Fig.|4]i. 

Note also that the use of y^ in Eq. ( ISSl l implies a slight ex- 
cess (~ (ys - ysWrlx^) with respect to the accurate numeri- 
cal results since ys decreases with the frequency (see Fig. |4|. 
In this representative test the excess is w 0.1 mK, but obvi- 
ously its value depends on the amplitude of free-free distortion 
(other than on the frequency). This error is clearl y negligible 
for the an alysis of current data (see e.g. [S alvaterra & Buriganal 
( l2002allbl) : IZannoni et all ( f2008l) : ISingaletal.i 6200^ ). It is likely 
not so relevant for the analysis of future measures at A ^ 1 cm 
with accuracy comparable to that proposed for DIMES (iKogutI 
ll996HBurigana & Salvaterrall2003h . On the contrary, it could be 
relevant for a very accurate analysis of future measures at with 
precision comparable to that proposed for FIRAS II / I $ 1 cm 
(Fixsen & Mather 2002; Burigana et al. 2004; Mathe3 l2009l) or 
thought in ideas for possible future long w avelength experiments 
from the Moon'^ 'llll (Burigana et al.ll2007l) . This calls for a com- 
plete frequency and thermal history dependent treatment of the 
free-free distortion in the accurate analysis of future data of ex- 
treme accuracy. 



6. Discussion and conclusion 

We have described the fundamental numerical approach and, in 
particular, the recent update to recent NAG versions of a nu- 
merical code, KYPRIX, specifically written for the solution of 
the Kompaneets equation in cosmological context, aimed to the 
very accurate computation of the CMB spectral distortions un- 
der quite general assumptions. The recipes and tests described in 
this work can be useful to implement accurate numerical codes 
for other scientific purposes using the same or similar numerical 
libraries or to verify the vaUdity of different codes aimed at the 
same or similar problems. 

Specifically, we have discussed the main subdivisions of the 
code and the most relevant aspects about technical specifications 
and code implementation. After a presentation of the equation 
formalism and of the boundary conditions added to the set of or- 
dinary differential equations derived from the original parabolic 
partial differential equation, we have given details on the adopted 
space (i.e. dimensionless frequency) grid, on the output results. 




10000 



1000 
A(cm) 



100 



2.727 



'u 2.726 



2.725 



: r 
; \ 
'■- \ 






: \ 






; \ 






'-- \ 






; v 
; \- 













2.724 

1000.00 100.00 10.00 



1.00 



0.10 



0.01 



A(cm) 



http://www.lnf.infn.it/conference/moon07/Program.html 



http://sci.esa.int/science-e/www/object/index.cfm7fobjec tid=40925 



Figure 4. Top panel: ys as derived from Eq. ( |54] |. Bottom panel: 
comparison between the numerical (long dashes) solution and 
the analytical approximation represented by Eq. dSSl) . See the 
text for further details. 



on the accuracy parameters, and on the used integration rou- 
tines. The introduction of the time dependence of the ratio be- 
tween electron and photon temperatures and of the radiative 
Compton scattering term, both introducing integral terms in the 
Kompaneets equation, has been addressed in the specific context 
of the recent NAG versions by discussing the solution adopted 
to solve the various related technical problems. 

Some of the tests we carried out to verify the reliability, ac- 
curacy, and performance of the code are presented. 

We have compared the results of the update version of the 
code with those obtained with the original one, reporting some 
representative cases, and we have found an excellent agreement. 

Some specific quantitative tests are here reported. They indi- 
cate a very good accuracy in the energy conservation: for appro- 
priate choices of the code accuracy parameters, the fractional in- 
jected energy is conserved within an accuracy better than 0.05%, 
or, in other words, possible energy conservation violations are 
negligible in practice for theoretical predictions and for com- 
parison with current and future data. The time behavior of the 
electron temperature is found to be in excellent agreement with 
the results obtained with the original code version, in spite of the 
different schemes adopted to update the evolving electron tem- 
perature. These are important verifications that probe that the 
cuiTent implementation of the code KYPRIX circumvents the 
problem represented by the lost possibility of internally adapt- 
ing the solver of partial differential equation to make it directly 
able to include integrals of the solution vector exactly at the cur- 
rent time step. Also, the setting of the accuracy parameters of the 
solver is now less flexible. Of course, these features imply a cer- 
tain increase of complexity in the code implementation, usage, 
and in the setting of code parameters. In spite of these difficul- 
ties, thanks to the better treatment of various computational steps 
involving integrals of functions over finite intervals, the new ver- 
sion of KYPRIX achieves a good accuracy even in the treatment 
of very small distortions. 

We have described also the introduction of the cosmological 
constant in the terms controlling the general expansion of the 
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Universe in agreement with the current concordance model, of 
the relevant chemical abundances, and on the ionization history, 
from recombination to cosmological reionization. 
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Figures. Final (i.e. at z = 1090) distorted spectra for differ- 
ent treatments of the recombination process. We assume as ini- 
tial condition a Comptonization spectrum at z = 10'* character- 
ized by Afr/e; - 10^^- Top panel: initial condition (dashes), final 
spectrum in the case of instantaneous recombination at z = 1090 
following the fully ionized phase (three dots-dashes), final spec- 
trum in the case of gradual, more realistic modelizations of the 
recombination process (the two different treatments give results 
indistinguishable in these plots). Bottom panel: difference be- 
tween the final solution found in the case of instantaneous re- 
combination and of gradual, more modelizations. See the text 
for further details. 

While an extensive discussion of various possible applica- 
tions of these new features of the code, object of future work, 
is not the scope of this paper, we report here few representative 
examples. 

In Fig. |5] we show the effect of different treatments of the 
computation of the hydrogen and helium ionization fractions 
during the recombination epoch in the case of a relatively late 
spectral distortion. Three different treatments of the recombina- 
tion history are exploited: an instantaneous recombination fol- 
lowing the fully ionized phase and two gradual, more realis- 
tic modelizations of the recombination process characterized by 
different evolutions of ionization fractions. In one of these two 
cases only the evolution of the electron ionization fraction, is 
taken from the code RECFAST while the H and He ionization 
fractions are computed solving the system of Saha equations; in 
the other case we directly use the complete output of the code 
RECFAST (see Sect.|43]and the caption of Fig.|5]for more de- 
tails). The results obtained in these two cases are identical in 
practice, while the assumption of full ionization up to z = 1090 
in the instantaneous recombination case produces a certain over- 
estimate (~ 10%) of the long wavelength excess due to (essen- 
tially free-free) photon production as well as a small amplifi- 
cation (~ 1%) of the relaxation of the initial Comptonization 
spectrum towards a Bose-Einstein like spectrum (obviously in- 
effecient, at the considered redshifts). With the adopted numer- 
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Figure 6. Distorted spectra predicted for the considered reion- 
ization models: initial blackbody spectrum (dots) at z = 20, 
almost corresponding to the beginning of the heating phase in 
the considered model, final (at z = 0) spectrum in the case of 
constant (f> (=147) (dashes) and for the realistic implementation 
of the ionization and thermal history in the suppression model 
(solid line). See the text for further details. 



ical accuracy (ACC=10 '2) and grid points (NPTS = 30001), 
the required CPU time on the used IBM platform is of about 40 
minutes, with less than 20% differences between the three dif- 
ferent cases. 

Fig- m reports the results derived for a particular astrophys- 
ical scenario of cosmological reionization. We adopt here the 
evolution of the electron ionization fraction Xe (the H and He 
ionization fractions being then computed solving the system 
of Saha equations) and of electro n temperature predicte d in 
the suppression model presented in [Schneider et al.l (l2008 f) and 
iBuripana et al.l (l2008h that implies a Thomson optical depth, 
T ^ 0.1017, almost in agreement with recent WMAP data (for 
consistency, we adopt in these computations the same cosmolog- 
ical parameters used in the quoted works). For comparison, we 
exploit a simple case with full ionization and a constant ratio, 0, 
between the electron temperature and the radiation temperature 
which is chosen to have the same amount of energy injected in 
the plasma {Aer/e, - 7.6 x 10"^) obtained in the above realistic 
model. As expected by construction, the two very different histo- 
ries produce the same high frequency (Comptonization like) dis- 
tortion but the long wavelength region is very different because 
of the contribution of free-free emission 'I The diff^er ence of this 
term in the two cases is mainly characterized by the evolution of 
the product (p^^^^xl ™d by the fact that free-free is more efficient 
at higher redshifts in the simple case than in the realistic one (see 
the caption of Fig.|6]for more details). With the adopted numer- 
ical accuracy (ACC=10"'^) and grid points {NPTS = 30001), 
the required CPU time on the used IBM platform is of few hours 



In both cases, the density contrast in the intergalactic medium 
associated to the formation of cosmic structures has been not included. 
The excess at low frequencies due to the free-free distortion should be 
then considered as lower limit in both cases, since the computation has 
been performed in the "averaged density" approximation. Therefore, a 
correction factor » <nj> / <ng>- > 1, coming from a proper inclusion 
of the treatment of density contrast in the intergalactic medium, should 
be apphed to the free-free term. 
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for the realistic model and about one hour in the simple model 
with a constant 0. 

Finally, we have focussed on some properties of the free-free 
distortions, relevant for the long wavelength region of the CMB 
spectrum, by checking that the new code version, as the original 
one, very accurately recovers the existing analytical approxima- 
tions in their limit of validity. We have also discussed the rele- 
vance of accurate computations able to improve the simple treat- 
ment based on the approximation with a frequency independent 
free-free distortion parameter 

All the tests done demonstrates the reliability and versatility 
of the new code version and its very good accuracy and appli- 
cability to the scientific analysis of current CMB spectrum data 
and of those, much more precise, that will be available in the 
future. 
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